Load libraries

### The packages (repetition) 
   library(here)
   library(sf)
   library(tmap)
   library(predicts)
   library(terra)
   library(viridis)
   #library(rgl) # not needed?
   #library(unmarked) # not needed?
   #library(spatstat) # not needed?

## check that the package raster is detached if you work with terra!
## However, it is later needed for the stack
  #library(raster) # not needed?
# detach(package:raster)

## only if you want to run MaxEnt from R ##

#install.packages('rJava') 
  library(rJava)
## please check if you get an error message
## If yes:
## 1. check if Java is installed on your computer
## if not: Download open jdk (https://jdk.java.net/23/),
## download the zip file and unpack to folder 'Java' under
## 'Programme' or 'program files'.
## Note: don’t use 'program files (86x)' - this may be the cause of error

## Then, set the system variable on your computer, -> Systemsteuerung, and on the search bar
## search for 'Umgebungsvariable', and then add a JAVA_HOME path to the path
## where Java is located on your Computer.


## Or set the directory of your Java location BEFORE loading the library:
## check with 

#Sys.getenv('JAVA_HOME') ## where the R-package is looking for the java.dll file
     ## set it to the correct path:
#Sys.setenv(JAVA_HOME='C:\\Program Files (x86)\\Java\\jre1.8.0_1447\\bin') 
     ## and recall the library
#library(rJava)

## http://www.r-bloggers.com/how-to-load-the-rjava-package-after-the-error-java_home-cannot-be-determined-from-the-registry/ 

Set the workspace

### The workspace  (repetition) 
getwd() # you can also use the package 'here()'

## [1] "C:/Users/wenzler/Documents/GitHub/d6_teaching_collection/R/Course3_a_presence_only_sdm"

root_wd <- here::here() 
# relative to work-wd

maps_wd <- paste(root_wd,"/","data/data_borneo/geo_raster_current_asc",sep='') # or:
maps_wd_fut <- paste(root_wd,"/","data/data_borneo/geo_raster_future_asc",sep='')
vecs_wd <- here::here("data","data_borneo","geo_vector") # shapefile
recs_wd <- here::here("data","data_borneo","animal_data") # location data

## the output folder should have been created by you during Tutorial 2 'R goes spatial'. 
## It should exist 
output_wd <- here::here("output")

# It should contain the hillshade.asc

setwd(output_wd) ## set to the OUTPUT folder!
getwd() # check

## [1] "C:/Users/wenzler/Documents/GitHub/d6_teaching_collection/output"

Load spatial data and define the CRS

### Data import - Load spatial data (repetition) 

ras1 <- terra::rast(x = paste(maps_wd,'/bio_asc_01.asc',sep=''))
# assign the projection (crs - coordinate reference system)
# ras1@crs <- CRS("+proj=longlat +datum=WGS84") ## not needed in terra!

ras24 <- terra::rast(here("data","data_borneo", "geo_raster_current_asc", "bio_asc_24.asc")) #DEM

ras42 <- terra::rast(x = here("data", "data_borneo", "geo_raster_current_asc", "bio_asc_42.asc")) # land use

#hillsh <- terra::rast(x = paste(maps_wd,'/borneo_hillshade.asc',sep='')) 

### Raster stack of environmental predictors (repetition) 

# the list.files command is very helpful to check what is in the folders
# use 'pattern' for searching for special file types, here .asc-files:
files <- list.files(path= maps_wd, pattern='.asc$',
                    full.names=TRUE )
files # these are just names! To load them as spatial objects, use rast()

## [1] "C:/Users/wenzler/Documents/GitHub/d6_teaching_collection/data/data_borneo/geo_raster_current_asc/bio_asc_01.asc"      
## [2] "C:/Users/wenzler/Documents/GitHub/d6_teaching_collection/data/data_borneo/geo_raster_current_asc/bio_asc_21.asc"      
## [3] "C:/Users/wenzler/Documents/GitHub/d6_teaching_collection/data/data_borneo/geo_raster_current_asc/bio_asc_22.asc"      
## [4] "C:/Users/wenzler/Documents/GitHub/d6_teaching_collection/data/data_borneo/geo_raster_current_asc/bio_asc_24.asc"      
## [5] "C:/Users/wenzler/Documents/GitHub/d6_teaching_collection/data/data_borneo/geo_raster_current_asc/bio_asc_27.asc"      
## [6] "C:/Users/wenzler/Documents/GitHub/d6_teaching_collection/data/data_borneo/geo_raster_current_asc/bio_asc_42.asc"      
## [7] "C:/Users/wenzler/Documents/GitHub/d6_teaching_collection/data/data_borneo/geo_raster_current_asc/borneo_hillshade.asc"

## note that I only load 4 rasters:
# predictors <- c(files[c(9,12,22,24)]) # for full data set
predictors <- terra::rast(x = c(files[c(1,5,6,7)]) ) # for github repository data

predictors

## class       : SpatRaster 
## dimensions  : 1425, 1423, 4  (nrow, ncol, nlyr)
## resolution  : 0.008333334, 0.008333334  (x, y)
## extent      : 108.3341, 120.1925, -4.374013, 7.500988  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84) 
## sources     : bio_asc_01.asc  
##               bio_asc_27.asc  
##               bio_asc_42.asc  
##               borneo_hillshade.asc  
## names       : bio_asc_01, bio_asc_27, bio_asc_42, borneo_hillshade

terra::plot(predictors, col = viridis::viridis(100)) ## might take some time depending on your computer

## what are the differences? Have a look also at the data description file in the data folder.

### Read in some Shapefiles (repetition) 
## package sf (new and maintained)

## Borneo outline polygon
Borneo_shp_sf <- sf::st_read(dsn = vecs_wd, 
                      layer = "borneo_admin",
                      stringsAsFactors = FALSE)[,c(1:3,5,7,17,18)]

## Reading layer `borneo_admin' from data source 
##   `C:\Users\wenzler\Documents\GitHub\d6_teaching_collection\data\data_borneo\geo_vector' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 158 features and 18 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 108.5986 ymin: -4.943054 xmax: 119.2697 ymax: 7.380556
## Geodetic CRS:  WGS 84

# Protected areas (PA) polygon
PA_shp_sf <-  sf::st_read(dsn = vecs_wd, 
                   layer = "Bor_PA",
                   stringsAsFactors = FALSE)[, c(1:4)]

## Reading layer `Bor_PA' from data source 
##   `C:\Users\wenzler\Documents\GitHub\d6_teaching_collection\data\data_borneo\geo_vector' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 255 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 108.5969 ymin: -4.18356 xmax: 119.6176 ymax: 9.911229
## Geodetic CRS:  WGS 84

# Rivers lines
River_shp_sf <- sf::st_read(dsn = vecs_wd, 
                     layer = "sn_100000",
                     stringsAsFactors = FALSE)

## Reading layer `sn_100000' from data source 
##   `C:\Users\wenzler\Documents\GitHub\d6_teaching_collection\data\data_borneo\geo_vector' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 396 features and 4 fields
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: 108.9454 ymin: -3.787917 xmax: 118.7879 ymax: 6.464583
## Geodetic CRS:  WGS 84

### The spatial point data (species records) (repetition) 

# filename
spec_pt_filename <- list.files(path = recs_wd, 
                               pattern = ".csv",
                               full.names = TRUE)#paste(recs_wd,'/','viverra_tangalunga.csv', sep='')
spec_pt_filename

## [1] "C:/Users/wenzler/Documents/GitHub/d6_teaching_collection/data/data_borneo/animal_data/DHOsim.csv"              
## [2] "C:/Users/wenzler/Documents/GitHub/d6_teaching_collection/data/data_borneo/animal_data/MyNewSpecies.csv"        
## [3] "C:/Users/wenzler/Documents/GitHub/d6_teaching_collection/data/data_borneo/animal_data/Pardofelis_marmorata.csv"
## [4] "C:/Users/wenzler/Documents/GitHub/d6_teaching_collection/data/data_borneo/animal_data/PPLsim.csv"              
## [5] "C:/Users/wenzler/Documents/GitHub/d6_teaching_collection/data/data_borneo/animal_data/RIVERsim.csv"            
## [6] "C:/Users/wenzler/Documents/GitHub/d6_teaching_collection/data/data_borneo/animal_data/viverra_tangalunga.csv"

# read the file
sp_recs <- read.csv(file = spec_pt_filename[1], # decide which species you want to use, here 1st in list
                    header=TRUE, sep=',')

#convert it to spatial object (sf here)
sp_recs_sf <- sf::st_as_sf(x = sp_recs, 
                       coords = c("long","lat"), # columns  for the coordinates
                       crs = 4326, # define crs, 4326 is the EPSG code
                       sf_column_name = "geometry",
                       remove=F) # sf needs a geometry column and you have to name it

Plot data overview

### Plot to get an impression 

plot(ras42, col=grey.colors(20))
plot(PA_shp_sf$geometry,border='green', lwd=1.8, add=T)
plot(sp_recs_sf$geometry, pch= '*',cex=1,col='deeppink',add=T)
text(112.3, 6, 'Starting with SDMs', cex=1.2, col= 'red')

Checking for multicollinearity

## workaround for slow computers. First, aggregate the 1 km² resolution into 50*50 km  
agg_pred <- terra::aggregate(x=predictors,fact=50,FUN=mean)
plot(agg_pred)

terra::pairs(agg_pred, method = 'spearman')

## Plot the differences (residuals) between the rasters:
diff <- terra::focalPairs(x = agg_pred, w = 3, 'pearson', na.rm = TRUE)
plot(diff)

# not_run on slow computers:
#pairs(predictors)

Sampling bias correction

We are now creating a very rudimentary bias correction file. Let’s assume that doing field work in the province of Sabah, Malaysia, is easy and the field sampling has concentrated on these areas. Further, it is easier to sample in lowlands and flat slopes. Usually, protected areas (PAs) are also well studied areas, where the sampling intensity is high. And let’s assume that it was not possible to get a research permit for Indonesia and that therefore very few field studies had been conducted here. Please note: this is a completely fictive example! You should have the information about the sampling bias.

That is, we create a file that is representing these differences in sampling effort. This is easiest to do with raster manipulation. We first create a raster from our shapefile with the country outlines.

head(Borneo_shp_sf)

## Simple feature collection with 6 features and 7 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 115.8098 ymin: 4.27081 xmax: 119.2697 ymax: 7.013332
## Geodetic CRS:  WGS 84
##   ID_0 ISO   NAME_0 NAME_1     NAME_2 Shape_Leng Shape_Area
## 1  158 MYS Malaysia  Sabah Lahad Datu  6.7475367 0.62106012
## 2  158 MYS Malaysia  Sabah      Papar  1.6727672 0.10065591
## 3  158 MYS Malaysia  Sabah  Penampang  0.9548857 0.03951545
## 4  158 MYS Malaysia  Sabah Pensiangan  4.0257859 0.49725547
## 5  158 MYS Malaysia  Sabah      Pitas  2.9223456 0.11955352
## 6  158 MYS Malaysia  Sabah      Ranau  2.4815704 0.24027591
##                         geometry
## 1 MULTIPOLYGON (((118.2299 4....
## 2 MULTIPOLYGON (((116.0144 5....
## 3 MULTIPOLYGON (((116.0477 5....
## 4 MULTIPOLYGON (((116.5357 5....
## 5 MULTIPOLYGON (((117.2762 6....
## 6 MULTIPOLYGON (((116.9182 6....

plot(Borneo_shp_sf[,3]) # column 3 is NAME_0 = main country names

We create a raster that gets the value of 1 for Malaysia (and Brunei) and a value of 0.01 (MaxEnt cannot deal with 0 -> see lecture) for the rest. We use one of our rasters as mask for cell size resolution and cropping extent:

### *Rasterize* the shapefiles
## command rasterize() 
## field = 1 means all raster cells get value = 1, other cell values are set to 0
## background sets all other cells to a chosen number


Mal_ras <- terra::rasterize(x=Borneo_shp_sf[Borneo_shp_sf$NAME_0 %in%
                       c("Brunei","Malaysia"),], y=ras24, field=1,
                       background=0.1)

plot(Mal_ras)

# almost same as not plotting Indonesia:
# NotInd_ras <- rasterize(x=Borneo_shp_sf[Borneo_shp_sf$NAME_0 !=
#                           "Indonesia",], y=ras24, field=1, background=0.01)

We do the same only for the province of Sabah and the Protected Areas

Sab_ras <- terra::rasterize(x=Borneo_shp_sf[Borneo_shp_sf$NAME_1 == 'Sabah',],
                     y=ras24, field=1, background=0.1)

PA_ras <- terra::rasterize(x=PA_shp_sf, y=ras24, field=1, background=0.1)


### Crosscheck with plot
par(mfrow=c(1,2))
plot(Sab_ras)
plot(PA_ras) #n.b.: only 1 and small values, no NA

par(mfrow=c(1,1))

And we do the same for all elevations below 300 m:

ras24_300 <- ras24 <= 300
plot(ras24_300)

Now we add them all up. This means, we rank the sampling intensity, i.e. sampling in protected areas in Sabah below 300 m gets the highest value

bias_tmp <- ras24_300 + Sab_ras + Mal_ras + PA_ras
bias_tmp

## class       : SpatRaster 
## dimensions  : 1425, 1423, 1  (nrow, ncol, nlyr)
## resolution  : 0.008333334, 0.008333334  (x, y)
## extent      : 108.3341, 120.1925, -4.374013, 7.500988  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84) 
## source(s)   : memory
## varname     : bio_asc_24 
## name        : bio_asc_24 
## min value   :        0.3 
## max value   :        4.0

plot(bias_tmp)

Let’s standardize between 0 and 1 (and then set the zeros to 0.01 again)

# get the maximum value in the layer, standardize it and round to two decimals
maxval <- max(values(bias_tmp),na.rm=T)
bias_tmp2 <- bias_tmp/maxval
bias_tmp3 <- round(bias_tmp2, digits=2)

## same as in one go:   
bias1 <- round(bias_tmp/max(values(bias_tmp),na.rm=T),digits=2)

table(values(bias1)) 

## 
##   0.08    0.3   0.33   0.53   0.55   0.75   0.78      1 
##  89415 107956 423878  39157 146727   5611  49196   4634

# table(bias1@data@values) # is doing the same

### in case of having 0 somewhere:
bias1[values(bias1 == 0)] <- 0.01 # because MaxEnt
                                  # does not take 0!
table(values(bias1))

## 
##   0.08    0.3   0.33   0.53   0.55   0.75   0.78      1 
##  89415 107956 423878  39157 146727   5611  49196   4634

Now compare the rasters - do they have the same extent, cell size etc?

terra::compareGeom(ras24,bias1) # the same?

## [1] TRUE

check for the length of the NODATA values:

# Are they really the same???
length(which(!is.na(values(bias1))))

## [1] 866574

length(which(!is.na(values(ras24))))

## [1] 866574

Nice! Now plot it:

### Plot the bias file
plot(bias1, col = viridis(7)) #or col = grey.colors(7)

## save into Output folder, which is your working directory: 
## check with getwd() if you're not sure of wd

## in case you want it aggregated
# bias_agg <- aggregate(bias1, fact=50,FUN=mean)
# writeRaster(bias_agg, filename="bias_agg.asc", datatype='ascii',
#            overwrite=TRUE,NAflag=-9999)

## since I used setwd() to the output folder, it is automatically stored there,
## if not, paste the path to the filename, i.e. 
## filename=paste0(output_wd,"/bias_2023.asc") #or # here(output_wd,"bias_2023.asc")
getwd()

## [1] "C:/Users/wenzler/Documents/GitHub/d6_teaching_collection/R/Course3_a_presence_only_sdm"

terra::writeRaster(bias1, filename="bias_2025.asc",
            overwrite=TRUE,NAflag=-9999)

dir.create(path = file.path(output_wd, "MaxEntRes"), showWarnings = FALSE)
dir.create(path = file.path(output_wd, "/MaxEntRes", "nobias_current"),showWarnings = FALSE)
dir.create(path = file.path(output_wd, "/MaxEntRes", "nobias_future"),showWarnings = FALSE)
dir.create(path = file.path(output_wd, "/MaxEntRes", "withbias_current"),showWarnings = FALSE)
dir.create(path = file.path(output_wd, "/MaxEntRes", "withbias_future"),showWarnings = FALSE)

select predictors

Here you can select which predictors you want to use. Please remember from the above exercise that current and future maps need to carry the same name, but in different folders, and that you must have the future maps for those variables that you included to fit the model to current conditions! That means, if you have fitted the model to current temperatures, you need to have maps of future temperatures as well (with the same name(s)) if you want to project the model to future climates.
The following example only builds on those files that are available in the repository, that is, the example runs on files bio_asc_01, bio_asc_21, bio_asc_22, bio_asc_24, bio_asc_27 and bio_asc_42. If you use the full dataset provided by the external link, please make sure that your predictor variables correspond.

files # files <- list.files(path= maps_wd, pattern='.asc$',full.names=TRUE )

## [1] "C:/Users/wenzler/Documents/GitHub/d6_teaching_collection/data/data_borneo/geo_raster_current_asc/bio_asc_01.asc"      
## [2] "C:/Users/wenzler/Documents/GitHub/d6_teaching_collection/data/data_borneo/geo_raster_current_asc/bio_asc_21.asc"      
## [3] "C:/Users/wenzler/Documents/GitHub/d6_teaching_collection/data/data_borneo/geo_raster_current_asc/bio_asc_22.asc"      
## [4] "C:/Users/wenzler/Documents/GitHub/d6_teaching_collection/data/data_borneo/geo_raster_current_asc/bio_asc_24.asc"      
## [5] "C:/Users/wenzler/Documents/GitHub/d6_teaching_collection/data/data_borneo/geo_raster_current_asc/bio_asc_27.asc"      
## [6] "C:/Users/wenzler/Documents/GitHub/d6_teaching_collection/data/data_borneo/geo_raster_current_asc/bio_asc_42.asc"      
## [7] "C:/Users/wenzler/Documents/GitHub/d6_teaching_collection/data/data_borneo/geo_raster_current_asc/borneo_hillshade.asc"

## Let's select the first 6 to eliminate the hillshade
predictors_selected <- terra::rast(x = c(files[c(1:6)]) ) 

Maxent without sampling bias correction

Here, the whole area that is given by the raster is sampled equally (background). Predictors are used as stack in the following MaxEnt command. Also, don’t forget to mark categorical variables and to not select variables which do not contain biological information, like the hillshade, as we did above.

xm_nobias_current <- predicts::MaxEnt(
                            x = predictors_selected, 
                            p = sp_recs[,2:3], # x and y coordinates
                            factors = 'bio_asc_42', # raster names which are factors
                            nbg = nrow(sp_recs) * 10, # number of background points, here ten times
                            path = file.path(output_wd, "/MaxEntRes", "nobias_current"), # output path
                            args = c("replicates=3", 
                                     "outputformat=logistic", 
                                     "responsecurves") # number of replicates and output format
                            )

saveRDS(xm_nobias_current, file.path(output_wd, "/MaxEntRes", "nobias_current", "maxent_nobias_current.rds"))

pred_no_bias_current <- terra::mean(predicts::predict(xm_nobias_current, predictors_selected))

terra::writeRaster(pred_no_bias_current, 
            file.path(output_wd, "/MaxEntRes", "nobias_current", "species_nobias_cur_avg.asc"),
            overwrite = TRUE)

## a little sneak preview:
terra::plot(pred_no_bias_current, col=viridis(100))

xm_nobias_future <- predicts::MaxEnt(
                            x = predictors_selected, # need to rereun (fit) MAxEnt to current conditions
                            p = sp_recs[,2:3], # x and y coordinates
                            factors = 'bio_asc_42', # raster names which are factors
                            nbg = nrow(sp_recs) * 10, # number of background points, ten times
                            path = file.path(output_wd, "/MaxEntRes", "nobias_future"), # output path
                            args = c("replicates=3", 
                                     "outputformat=logistic", 
                                     "responsecurves",
                            paste0("projectionlayers=", maps_wd_fut)) 
                            )

saveRDS(xm_nobias_future, file.path(output_wd, "/MaxEntRes", "nobias_future", "maxent_nobias_future.rds"))

## sneak preview - ascii file should be in folder:
## check for the average '_avg' file out of all single repetitions:
fut_avg_nobias <- list.files(here::here("output","MaxEntRes","nobias_future"),pattern='_avg.asc$',
                    full.names=TRUE )
fut_avg_nobias

## [1] "C:/Users/wenzler/Documents/GitHub/d6_teaching_collection/output/MaxEntRes/nobias_future/species_geo_raster_future_asc_avg.asc"

sneak <- terra::rast(x = fut_avg_nobias)
terra::plot(sneak,col= viridis(100))

## NB: the future plot might not look very different from the current one,
## because not many changing predictor variables were used.

## create a residual plot subtracting current from future
## as the output of MaxEnt is a RasterLAyer, 
## it needs to be transformed to a terra::SpatRaster 
res_plot <- sneak - pred_no_bias_current
plot(res_plot, col= inferno(100))

MaxEnt with sampling bias correction

## The no. of bg points was increased to * 50 to avoid conversion issues
## should be done consistently also across nobias-runs, but not done
## here to save computation time.
random_pts <- terra::spatSample(x = bias1, 
                                size = nrow(sp_recs) * 50, 
                                method = "weights", 
                                na.rm=TRUE, 
                                as.points=TRUE)

xm_withbias_current <- predicts::MaxEnt(
                            x = predictors_selected, 
                            p = sp_recs[,2:3], # x and y coordinates
                            factors = 'bio_asc_42', # raster names which are factors
                            a = random_pts, # bias - background points
                            path = file.path(output_wd, "/MaxEntRes", "withbias_current"), # output path
                            args = c("replicates=3", 
                                     "outputformat=logistic", 
                                     "responsecurves") # number of replicates and outputformat
                            )

saveRDS(xm_withbias_current, file.path(output_wd, "/MaxEntRes", "withbias_current", "maxent_withbias_current.rds"))

pred_with_bias_current <- terra::mean(predicts::predict(xm_withbias_current, predictors_selected))

terra::writeRaster(pred_with_bias_current, 
            file.path(output_wd, "/MaxEntRes", "withbias_current", "species_bias_cur_avg.asc"),
            overwrite = TRUE)

## a little sneak preview:
terra::plot(pred_with_bias_current, col=viridis(100))

xm_withbias_future <- predicts::MaxEnt(
                            x = predictors_selected, 
                            p = sp_recs[,2:3], # x and y coordinates
                            factors = 'bio_asc_42', # raster names which are factors
                            a = random_pts, # bias - background points
                            path = file.path(output_wd, "/MaxEntRes", "withbias_future"), # output path
                            args = c("replicates=3", 
                                     "outputformat=logistic", 
                                     "responsecurves",
                            paste0("projectionlayers=", maps_wd_fut))
                            )

saveRDS(xm_withbias_future, file.path(output_wd, "/MaxEntRes", "withbias_future", "maxent_withbias_future.rds"))

Accessing MaxEnt results

Have a look at the MaxEnt layers and results produced in the folder MaxEntRes. MaxEnt saves the average of x repetitions (remember, we ran 3 repetitions) in a file with ’_avg.asc’ at the end of the name.Please take care that you have results in each of the four folders under MaxEntRes, and that the folders are built in the way we did it above.

# new argument 'recursive' means, that also subfolders are checked!
infiles <- list.files(path=paste(output_wd,'/MaxEntRes',
                                 sep=''),pattern='_avg.asc$',
                      full.names=TRUE,recursive=TRUE )
infiles

## [1] "C:/Users/wenzler/Documents/GitHub/d6_teaching_collection/output/MaxEntRes/nobias_current/species_nobias_cur_avg.asc"            
## [2] "C:/Users/wenzler/Documents/GitHub/d6_teaching_collection/output/MaxEntRes/nobias_future/species_geo_raster_future_asc_avg.asc"  
## [3] "C:/Users/wenzler/Documents/GitHub/d6_teaching_collection/output/MaxEntRes/withbias_current/species_bias_cur_avg.asc"            
## [4] "C:/Users/wenzler/Documents/GitHub/d6_teaching_collection/output/MaxEntRes/withbias_future/species_geo_raster_future_asc_avg.asc"

# example
#[1] "c:/Users/kramer/Dropbox/.../MaxEntRes
#     /nobias/river_dummy_avg.asc"
#[2] "c:/Users/kramer/Dropbox/.../MaxEntRes
#     /nobias/river_dummy_FutureMaps_avg.asc"
#[3] "c:/Users/kramer/Dropbox/.../MaxEntRes
#     /withbias/river_dummy_avg.asc"
#[4] "c:/Users/kramer/Dropbox/.../MaxEntRes
#     /withbias/river_dummy_FutureMaps_avg.asc"

Stack results file and make usual plot

## IF you have a slow computer, 
## please use the workaround with function aggregate() from above:
# me_stack <- aggregate(c(infiles[c(1:4)]),fact= 50, FUN = mean )

me_stack <- c(rast(infiles[[1]]),
              rast(infiles[[2]]),
              rast(infiles[[3]]),
              rast(infiles[[4]]))

# name sequence according to infiles list above
names(me_stack) <- c('curr_noBias','fut_noBias',
                     'curr_withBias','fut_withBias')

#par(mar=c(4,4,3,3), oma=c(2,2,2,2))

plot(me_stack,col=viridis(100)) # note: I might use a different example than you did in MaxEnt

Check the range of the suitability values

boxplot(me_stack,layers = c(1,3,2,4),notch=T,outline=F)

How would you interpret the results?

Threshold selection

You should have results for the runs with a bias file and without in a .csv file. When you run MaxEnt via the GUI, you receive only two results files, because the model is fitted to the data of the current condition; in R you will receive four files, simply because running MaxEnt with R for future maps is a re-fit of the model to current conditions. The files should be identical, so it is enough to only save the current results file.

# Use maxentResults.csv to access threshold value
# column: 10th percentile training presence logistic threshold

me_res <- list.files(path=paste(output_wd,'/MaxEntRes',sep=''),
                     pattern='maxentResults.csv',
                     full.names=TRUE,recursive=TRUE )
me_res

## [1] "C:/Users/wenzler/Documents/GitHub/d6_teaching_collection/output/MaxEntRes/nobias_current/maxentResults.csv"  
## [2] "C:/Users/wenzler/Documents/GitHub/d6_teaching_collection/output/MaxEntRes/nobias_future/maxentResults.csv"   
## [3] "C:/Users/wenzler/Documents/GitHub/d6_teaching_collection/output/MaxEntRes/withbias_current/maxentResults.csv"
## [4] "C:/Users/wenzler/Documents/GitHub/d6_teaching_collection/output/MaxEntRes/withbias_future/maxentResults.csv"

# example
#[1] "c:/Users/kramer/Dropbox/.../MaxEntRes
#     /nobias/maxentResults.csv"
#[2] "c:/Users/kramer/Dropbox/.../MaxEntRes
#     /withbias/maxentResults.csv"

Now we do some elegant trick: we store both files as two objects in a list! Similar to stacking rasters, a bit….

## Please choose here:
store_res <- lapply(me_res,FUN=read.csv) ## store as list obj. when run from GUI
store_res <- lapply(me_res[c(1,3)],FUN=read.csv) ## when run from R

## Check the following out:
# head(store_res) # list object with 2 slots
# str(store_res)
# store_res[[1]]

We are searching for the column with the threshold.

## Extract the values at column '10th percentile' at position nrep +1
## (= nrow command) to access avg value of all repetitions at the lowest row

zename<-'X10.percentile.training.presence.Logistic.threshold' 
## Note! :
## In MaxEnt 3.4.1 the column name is:
## 'X10.percentile.training.presence.Logistic.threshold'
## with capital L in 'logistic', in older versions with 'l'.
zename

## [1] "X10.percentile.training.presence.Logistic.threshold"

Now check which column number that is….

zecol <-which(colnames(store_res[[1]]) == zename) 
zecol ## column 51 in my case

## [1] 51

… and extract the threshold value :

## a little head twister:
## access each element of your list, which is a whole dataset with store_res[[1]]
## and inside this list object = data.frame, access the element of
## the last row (nrow), and the column zecol containing the threshold value,
## because the results are stored for each repetition, and the last column contains
## the average value
t_noBias   <- store_res[[1]][nrow(store_res[[1]]),zecol]
t_withBias <- store_res[[2]][nrow(store_res[[2]]),zecol]

t_noBias #the thresholds

## [1] 0.3625

t_withBias

## [1] 0.387

Now apply the thresholds to your MaxEnt output maps of relative probability. For that, store the thresholds in a vector, and double them for bias/ nobias, because the threshold stays the same for current and future maps (see above):

all_bias <- rep(c(t_noBias,t_withBias), each=2)
all_bias

## [1] 0.3625 0.3625 0.3870 0.3870

Now set the thresholds to the stack. If you don’t remember the stack (me_stack), plot it again. Test different thresholds and see how the predictions change:

binary_thresh <- me_stack >= all_bias
plot(binary_thresh)

test05 <- me_stack >= 0.5 # check result with fixed threshold of default cut-off
plot(test05)

#testSeq <- me_stack >= seq(0.1, 0.7, length = 4)
#plot(testSeq)

Analyses

Now let’s check how global change affects the protected areas. Does the suitability of these areas increase or decrease in the future? Extract relative probability values (habitat suitability) inside protected areas PAs and make a boxplot for comparison.

## check if overlay correct - only works when you have NOT aggregated the maps
terra::compareGeom(PA_ras,me_stack[[1]]) 

## [1] TRUE

## extract
ex <- which(values(PA_ras)==1) ##remember PA_ras == 1 where the PAs are
## ex # gives Pointer to elements/ Index
ex_stack <- terra::extract(x=me_stack, y=ex)
head(ex_stack)

##    curr_noBias  fut_noBias curr_withBias fut_withBias
## 1           NA          NA            NA           NA
## 2 0.0011197921 0.001269990  2.361548e-04  3.07386e-05
## 3 0.0004946960 0.000529586  1.733264e-04  4.35029e-05
## 4 0.0006508699 0.000678241  1.230861e-04  5.65338e-05
## 5 0.0003685336 0.000376808  1.551008e-04  3.38545e-05
## 6 0.0001942456 0.000186320  5.513565e-05  2.72841e-05

Plot the difference of habitat suitability in scenarios and PAs

boxplot(ex_stack,na.rm=T)

Calculate no. of cells (=area) of suitable habitat in PAs, that we defined via the binary threshold maps.

ex_binary <- terra::extract(x=binary_thresh, y=ex)
head(ex_binary)

##   curr_noBias fut_noBias curr_withBias fut_withBias
## 1          NA         NA            NA           NA
## 2       FALSE      FALSE         FALSE        FALSE
## 3       FALSE      FALSE         FALSE        FALSE
## 4       FALSE      FALSE         FALSE        FALSE
## 5       FALSE      FALSE         FALSE        FALSE
## 6       FALSE      FALSE         FALSE        FALSE

## In the following, we are summing the cells with 'TRUE'
colSums(ex_binary,na.rm=T) 

##   curr_noBias    fut_noBias curr_withBias  fut_withBias 
##         46892         53211         47519         33686

Interestingly: when correcting for sampling bias, the future scenarios have the lowest number of suitable habitat inside PAs (in my species):

curr_noBias fut_noBias curr_withBias fut_withBias
43285 63905 43430 39989

old not run

run from here again

## where does the animal corridor intersect rivers?
## in sf package
River_inter1 <- sf::st_intersection(x = River_shp_sf, y = sPath_sf) 

plot(sf::st_geometry(sPath_sf),lwd=3) # we need st_geometry to have a good plot
plot(PA_shp_sf[,1],border='grey',lwd=1,add=T)
plot(River_shp_sf[,1],col='blue',lwd=2,add=T)
plot(River_inter1[,1],col='red',lwd=4,add=T)
box()

## Which connections intersect the PAs?

## as PA_shp_sf has an invalid topology, we need the workaround with the old 
## {rgeos}package; https://github.com/r-spatial/sf/issues/1902
PA_inter1 <- sf::st_intersection(sf::st_make_valid(PA_shp_sf), sPath_sf) # old way

plot(sf::st_geometry(sPath_sf),lwd=3)  # we need st_geometry to have a good plot
plot(PA_shp_sf,border='blue', col='transparent', lwd=2,add=T)
lines(PA_inter1,col='red',lwd=4); box()

### Save new line as shapefile
sf::st_write(obj = sPath_sf, dsn = paste0(output_wd,'/costpath_line.shp'), delete_layer = TRUE)

## Deleting layer `costpath_line' using driver `ESRI Shapefile'
## Writing layer `costpath_line' to data source 
##   `C:/Users/wenzler/Documents/GitHub/d6_teaching_collection/output/costpath_line.shp' using driver `ESRI Shapefile'
## Writing 1 features with 3 fields and geometry type Line String.

################  end of the course #######################